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Abstract. We use 3-dimensional numerical simulations of self-gravitating compress- 
ible turbulent gas in combination with Lagrangian tracer particles to investigate the 
mixing process of molecular hydrogen (H2) in interstellar clouds. Tracer particles are 
used to represent shock-compressed dense gas, which is associated with H2. We deposit 
tracer particles in regions of density contrast in excess of 10jOo> where po denotes the 
mean density. Following their trajectories and using probability distribution functions, 
we find an upper limit for the mixing timescale of H2, which is of order 0.3 Myr. This 
is significantly smaller than the lifetime of molecular clouds, which demonstrates the 
importance of the turbulent mixing of H2 as a preliminary stage to star formation. 
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1. INTRODUCTION 

Turbulent mixing of chemical species in the interstellar medium (ISM) has profound 
consequences on the morphology, chemistry, cooling and shielding of molecular clouds 
(MCs). MCs are known to be the formation site of young stars. It is therefore 
a prerequisite for a successful theory of star formation to explain the physical and 
chemical environment in which stars are born. Measured velocity dispersions suggest 
that interstellar clouds exhibit internal supersonic compressible random motions, that 
are associated with turbulence (e.g., Larson 1981, Heyer & Brunt 2004). Kinetic energy 
power spectra have been measured (e.g., Padoan et al. 2006), which appear to be slightly 
steeper than the Kolmogorov (1941) spectrum of incompressible turbulence, even if 
intermittency corrections are considered (e.g., Kolmogorov 1962, She & Leveque 1994, 
Boldyrev 2002). Excellent summaries of empirical and theoretical aspects of interstellar 
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turbulence are presented in the review articles by Elmegreen & Scalo (2004), Scalo & 
Elmegreen (2004) and Mac Low & Klessen (2004). 

Compressions created by the formation of strong shocks facilitate the formation of 
molecular hydrogen, H2 (e.g., Hollenbach et al. 1971). Observations of MCs typically 
detect carbon monoxide as a tracer for H2, which is found to be ubiquitous throughout 
MCs and not merely in regions of high density. Turbulent transport of H2 from high 
density to low density gas may be an explanation for its relative homogeneity. 

A recent numerical investigation of interstellar turbulence by Glover & Mac 
Low (2007) confirms that turbulent compressions are a key mechanism for the rapid 
(1 — 2Myr) formation of H2. Their results further suggest that a significant fraction of 
H2 is located in low density gas, in quantities greater than can be formed directly in 
situ. The authors propose that H2 forms in shock-compressed sheets and filaments, and 
is subsequently transported to lower density regions. In order to test this hypothesis, 
we have performed simulations of interstellar turbulence that employ Lagrangian tracer 
particles to directly follow the trajectories of dense gas parcels. Tracer particles are 
deposited in high density regions, formed in a driven compressible turbulent medium. 
We then use statistical analysis to show that even under self-gravitating conditions, 
significant mixing of H2 occurs on short timescales (< 0.3 Myr). 

In section 2, we describe our numerical schemes and simulation setup. Section 3 
presents our results and discussion, and in section 4, our conclusions are summarized. 



2. NUMERICAL METHODS AND SIMULATION SETUP 

We solve the equations of compressible hydrodynamics including self-gravity on a 3- 
dimensional static grid of 256^ zones using the piecewise parabohc method (PPM) 
(Colella & Woodward 1984), implemented in the astrophysical code ENZO (O'Shea 
et al. 2004) with periodic boundary conditions. Density p, velocity v, total energy 
density pe, pressure P and gravitational potential $ are related through the equations 

|^ + V-(pv) = (1) 
dv 1 

— + (v • V)v = — VP - V$ + f (2) 
at p 

^(pe) + V ■ [v(pe + P)] = -pv ■ (V$) + p v f (3) 
A$ = AttGp , (4) 
where G is the gravitational constant. An isothermal equation of state, 

P = p(7 — l)u with u = e — -v^ , (5) 

approximated by 7 = 1.01 is used to close the equations of hydrodynamics. This is 
a crude, but reasonable approximation on a wide range of scales in both length and 
density in interstellar clouds. 
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We start from a uniform distribution of gas, initially at rest. In order to excite 
turbulent motions, we utilize a stochastic forcing term f, appearing as source term in 
equations (2) and (3) that supplies kinetic energy on the largest scales. The forcing term 
has a parabolic Fourier spectrum with k/k^ = [0.5, 1.5] centered on ko = 27r/L, where L 
is half of the size of the computational domain, and is evolved by an Ornstein-Uhlenbeck 
process (e.g., Eswaran & Pope 1988, Schmidt et al. 2006). In this study, we adjust our 
forcing amplitude such that the fluid reaches a representative rms Mach number of 3. 
We emphasize that our forcing is constructed such that we can regulate the relative 
strength of compressive modes (V x f = 0) with respect to solenoidal modes (V • f = 0). 
The ratio of kinetic energies at the injection scale, 

^ = :bttI — 

-"^sol ~ -"-^comp 

is approximately 0.8, which represents mostly solenoidal forcing. Note that x cannot be 
exactly 1, despite the utilization of purely solenoidal forcing, because supersonic flows 
always contain compressive modes. The consequences of varying the forcing amplitude 
and mode are discussed in Schmidt & Federrath (2007). 

Timescales are measured in units of the autocorrelation timescale T of the forcing, 
which was set equal to the turbulent crossing time on large scales for a Mach 3 turbulent 
medium. This means that after approximately 1 T, the gas reaches an rms Mach number 
Mrms ~ 3, which is maintained by the forcing for all times t > IT. Within < t < IT, 
self-gravity is kept deactivated, allowing the gas to reach a state of fully developed 
compressible turbulence. After that, at t = IT, Lagrangian tracer particles are placed 
at the centers of grid cells with density p/po > 10 (po is the mean density), representing 
highly compressed gas parcels. From this point on, we distinguish three routes of further 
evolving the simulation: 

1) no changes, self-gravity still deactivated (pure forcing) 

2) self-gravity activated, representing 64 Jeans masses (forcing + SG4) 

3) self-gravity activated, representing 4096 Jeans masses (forcing -|- SG64) . 

Case 1 (pure forcing) serves as a control run to compare to the self-gravitating cases 
2 and 3. The difference between 2 (SG4) and 3 (SG64) lies in the mean density po of 
the fluid, which makes up a domain that contains a total of 64 Jeans masses for case 2 
(SG4) and a total of 4096 Jeans masses for case 3 (SG64). The Jeans mass is defined 
by 

Mj = poA/ = Po^ , (7) 



\Gpo 

where Cg = a/ 7-P/ p is the sound speed. This means that the self-gravitating cases 
represent gravitationally unstable gas, and both cases would be subject to gravitational 
collapse in the absence of turbulent fluctuations and stochastic forcing. 

The Langrangian tracer particles are evolved every timestep (using a simple Euler- 
step) according to the Eulerian velocity, which is interpolated from the underlying grid 
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Figure 1. Spatial distribution of tracer particles at the instance of their deposition 
t = LOT in regions with p/po > 10. They represent shock-compressed gas, which is 
associated with H2. 



by means of a second order accurate (triangular shaped cloud) method at the current 
position of each tracer particle. We have also tried a first order (cloud in cell) and a 
third order (tricubic, Lekien & Marsden 2005) spatial interpolation, and a predictor- 
corrector-step for time integration, which yielded no statistically significant differences. 

3. RESULTS AND DISCUSSION 

Figure 1 shows the spatial distribution of tracer particles at t = LOT", the instant of 
their deposition, which serves as an initial condition for three different routes of further 
development (see section 2). Approximately 24000 tracer particles have been placed in 
regions of shock compressed gas with density contrast in excess of 10, representing H2 
in the following. In figure 2, we compare the evolution of the pure forcing run (left 
panel) to the 'strong' self-gravitating case SG64 (right panel) . The case of 'weak' self- 
gravity (SG4) is not shown here, because it is almost identical to the evolution of the 
pure forcing run. This is because the forcing still dominates the dynamics in this case, 
whereas self- gravity becomes the dominant force in case SG64. Visual inspection of 
SG64 in comparison to pure forcing reveals that mixing of tracer particles occurs for 
both cases, but mixing is less efScient in the self-gravitating case. This is to be expected, 
since self-gravity acts to confine dense cores, within which the turbulent kinetic energy 
is less than the potential energy. Tracer particles located in such regions are trapped 
by self-gravity. 

In order to quantify this behavior, we compare probability distribution functions 
(PDFs) of the Eulerian gas density (density defined on the grid) and PDFs of Lagrangian 
gas density (interpolated from the grid and monitored at the position of each tracer 
particle). In fully developed compressible turbulent flows, the PDF of density is close to 
a log-normal distribution (e.g., Padoan et al. 1997, Klessen 2000, Ostriker et al. 2001, Li 
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pure forcing forcing + self gravity 




Figure 2. Temporal evolution of the spatial distribution of Lagrangian tracer particles 
(from top to bottom). Left: simulation without self-gravity (pure forcing). Right: 
simulation including self-gravity (forcing + SG64). Self-gravity can prevent some of 
the tracer particles from mixing with the background gas. There are no snapshots for 
t > 1.2T for case SG64, because SG64 violates the resolution criterion for simulations 
including self-gravity (Truelove et al. 1997) shortly after t = 1.2 T, which only allows us 
to draw reliable conclusions up to this time. Case SG4 (not shown here) is statistically 
similar to the case of pure forcing, as discussed in the text. 



et al. 2004, Kritsuk et al. 2007) and given by 



Pv,m{s)ds 



27r(7 



exp 



2^^ 



ds 



(8) 
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Figure 3. Time sequence of ttie mass-weiglited density PDFs for the case without 
self-gravity (pure forcing). The Lagrangian PDF of the tracer particles (solid line) 
gradually approaches the Eulcrian PDF (dashed line). The dotted lines are least- 
squares- fits using the log-normal distribution given by equation (8). 

with s = ln(p/po), mean Sy^m and variance a^. The subscripts v and m denote volume- 
and mass-weighted distributions, respectively. Volume- and mass-weighted PDFs are 
related by = —Sm = — (e.g., Li et al. 2003). We make use of this fact to 
transform the volume-weighted log-normal Eulerian PDF into the corresponding mass- 
weighted PDF. This allows us to directly compare the Eulerian PDF to the Lagrangian 
PDF of the tracer particles, which is naturally a mass-weighted distribution. 

Figure 3 shows a time sequence of the mass- weighted PDFs of both Eulerian (dashed 
lines) and Lagrangian (solid lines) gas density for the case of pure forcing. At t = 1.00 T, 
the tracer particle distribution covers densities in excess of s ~ 2.3, which corresponds to 
logio[p/Po] ~ 1, as expected from their deposition criterion. The Eulerian distribution is 
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Figure 4. Same as figure 3, but for the self-gravitating case SG64. The Lagrangian 
PDF of the tracer particles (solid line), as well as the Eulerian PDF (dashed line) 
deviate from a log-normal distribution at high densities and never fall on top of each 
other. In contrast to the case of pure forcing, complete mixing is impossible due 
to gravitational confinement. However, also in this case, turbulent transport to low 
density regions is very efficient. The dotted lines are least-squares-fits using the log- 
normal distribution given by equation (8). 



in good agreement with a log-normal distribution, indicated by a fit using equation (8) 
(dotted line), which justifies our transformation from volume- to mass- weighted PDFs. 
Following the time sequence, it is evident that the tracer particle PDF converges to the 
Eulerian PDF. For t > 1.30 T, both Eulerian and Lagrangian distributions are (except 
for small fiuctuations in the wings of the distributions) statistically identical. This 
behavior indicates mixing of tracer particles to rarefied lower density regions. 

In contrast to the case of pure forcing, the PDFs from SG64 shown in figure 4 never 
reach a state where they resemble each other. Especially the high density wings, which 
are mostly affected by self-gravity deviate from a log-normal distribution. The Eulerian 
PDF tends to develop a power law tail for high density, in agreement with Klessen (2000). 
Although significant mixing into low density regions occurs within 0.1 T of the tracer 
particle deposition, some tracers never completely mix with the background gas, owing 
to the self-gravity of the gas that holds them within deep gravitational potential wells 
in high density regions. The onset of gravitational collapse manifests in the increasing 
probability for tracer particles to be found at very high density, which is indicated by 
the fiattening of the Lagrangian PDF towards high density. 
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Figure 5. Cumulative probability distribution of the density, which is monitored by 
Lagrangian tracer particles at t = 1.1 T for the non-sclf-gravitating case (pure forcing, 
solid line), the 'weak' self-gravitating case (SG4, dotted line) and the 'strong' self- 
gravitating case (SG64, dash-dotted line). The vertical dashed line marks the density 
p = lOpo, which has been the deposition threshold for tracer particles at t = 1.0 T. 
All tracer particles, which contribute to the left of this line have been mixed to lower 
density regions. Pure forcing and SG4 are almost identical, because turbulence still 
dominates the dynamics. For SG64, 70% of the tracer particles have been mixed to 
low density regions. For pure forcing and SG4, it is 95%. 



The fraction of tracer particles that have been mixed to regions of density smaller 
than their deposition density (p < lOpo) can be estimated using cumulative distribution 
functions (CDFs) of the Lagrangian density distribution. We compare CDFs of all 
three runs (pure forcing, SG4 and SG64) at time t = 1.1 T in figure 5. For the 
simulation without self-gravity and for SG4, approximately 95% of the tracer particles 
have been mixed to regions of density smaller than lOpo (vertical dashed line). In 
the case of 'strong' self-gravity (SG64), approximately 70% have been transported to 
low density. This quantifies the influence of self-gravity in run SG64 on the mixing 
efficiency compared to pure forcing and SG4 cases. Although a significant fraction of 
tracer particles remains at high density, most of the tracer particles have already been 
mixed into lower density regions within 0.1 T after their deposition, even in case SG64. 



Turbulent mixing of shock-compressed dense gas, which is associated with molecular 
hydrogen (H2) and followed by means of Lagrangian tracer particles has been 
investigated in self-consistent 3-dimensional numerical simulations of driven supersonic 
turbulence. Even for self-gravitating media, turbulent transport from high density cores 
to low density rarefied regions is very efficient. Using Lagrangian tracer particles allowed 
us to directly confirm and quantify the importance of turbulent transport as the main 



4. CONCLUSIONS 
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mechanism for H2 mixing in the intersteUar medium (see Glover & Mac Low 2007). We 
estimate that approximately 70% of the H2 formed in shock-compressed clumps and 
filaments can be mixed into lower-density gas within one tenth of a turbulent crossing 
time for a solenoidally driven turbulent medium. This represents an upper limit for 
the mixing timescale of 0.3 Myr for typical interstellar cloud conditions in the cold 
neutral medium (rms Mach number 3, cloud diameter lOparsec, sound speed 0.6kms~^, 
dimensions similar to Glover & Mac Low (2007)). However, interstellar turbulence is 
likely to be driven compressively (e.g., Federrath et al. 2007, and references therein), 
a situation that is markedly different from solenoidally driven turbulence (Schmidt & 
Federrath 2007) . This needs to be quantified in a follow-up study. 
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